Distinguishable Dyads - Cross-Sectional APIM

Distinguishable Dyads - Data: Simulated Dyads

print_df(head(df))

Distinguishable Dyads - Cross-Sectional APIM

(e.g., Kenny et al., 2006; Kenny & Cook, 1999; Kenny & Kashy, 2011)

Distinguishable Dyads - Preparing Data

df_apim <- df %>%
  # Reshaping to Actor-Partner format (4-field)
  reshape_dyadic_data(
    person_id = 'userID',
    dyad_id = 'coupleID',
    vars_to_reshape = 'communication'
  )  %>%
  mutate(
    # Optional: grand-mean centering
    communication_actor_gmc = communication_actor - mean(communication_actor),
    communication_partner_gmc = communication_partner - mean(communication_partner),
    
    # Create Dummy-Variables for male and female
    is_female = ifelse(gender == 1, 1, 0),
    is_male = ifelse(gender == 2, 1, 0),
  ) 

print_df(head(df_apim))

Distinguishable Dyads - Preparing Data

Distinguishable Dyads - Fitting the Model in BRMS

Model non-independence either via a shared random intercept or correlated residuals (not both). Posterior predictive checks, residual diagnostics and leave one out cross validation (e.g., via loo_compare()) may help inform which one to use.

For other families like ordinal regression with brms, one might work better than the other.

Distinguishable Dyads - Fitting the Model in BRMS

# Option A: Random-intercept model for non-independence

formula <- bf(
  satisfaction ~ 
    
    # Remove global intercept, introduce male and female intercepts.
    0 + is_male + is_female +
    
    # Actor effect for male and female
    communication_actor_gmc:is_male + 
    communication_actor_gmc:is_female + 
    
    # Partner effect for male and female
    communication_partner_gmc:is_male + 
    communication_partner_gmc:is_female +
    
    # Option 0: Random intercept for male and female (correlated)
    # NOT identifiable in the cross-sectional case:
    # (0 + is_male + is_female | coupleID)
    
    # Option A: Shared dyad effect - interpretable dyad-level random intercept. 
    (1 | coupleID), 
  
  # Two distinct residual variances for males vs females:
  sigma = ~ 0 + is_male + is_female
)

priors <- c(
  # prior(normal(2, 3), class = "Intercept"),
  prior(normal(0, 5), class = "b"),
  prior(student_t(3, 0, 1.5), class = "b", dpar = "sigma"),
  prior(student_t(3, 0, 2.5), class = "sd")
)

model_dist_apim_a <- brm(
  formula = formula, 
  data = df_apim,
  family = gaussian(link = identity),
  prior = priors,
  chains = 4,
  cores = 4,
  iter = 2000,
  warmup = 1000, 
  seed = 123,
  file = file.path('brms_cache', 'example1_dist_apim_a') # Cache the model
)

Distinguishable Dyads - Fitting the Model in BRMS

# Option B: Residual CS for non-independence (no random intercept)

formula <- bf(
  satisfaction ~ 
    
    # Remove global intercept, introduce male and female intercepts.
    0 + is_male + is_female +
    
    # Actor effect for male and female
    communication_actor_gmc:is_male + 
    communication_actor_gmc:is_female + 
    
    # Partner effect for male and female
    communication_partner_gmc:is_male + 
    communication_partner_gmc:is_female +
  
    # Option B: Model non-independence by using compound symmetry of the residuals. 
    # Aligns with the SEM conceptualization / literature
    cosy(gr = coupleID),
  
  # Two distinct residual variances for males vs females:
  sigma = ~ 0 + is_male + is_female
)

priors <- c(
  # prior(normal(2, 3), class = "Intercept"),
  prior(normal(0, 5), class = "b"),
  prior(student_t(3, 0, 1.5), class = "b", dpar = "sigma"),
  prior(beta(1, 3), class = "cosy")
)

model_dist_apim_b <- brm(
  formula = formula, 
  data = df_apim,
  family = gaussian(link = identity),
  prior = priors,
  chains = 4,
  cores = 4,
  iter = 2000,
  warmup = 1000, 
  seed = 123,
  file = file.path('brms_cache', 'example1_dist_apim_b') # Cache the model
)

Distinguishable Dyads - Check Model Convergence and Fit

Check Rhats and Effective Sample Sizes (ESS_tail and ESS_bulk) directly from the brms summary. Additionally you can (using Model B as an example):

rstan::check_hmc_diagnostics(model_dist_apim_b$fit)

Divergences:

0 of 4000 iterations ended with a divergence.

Tree depth:

0 of 4000 iterations saturated the maximum tree depth of 10.

Energy:

E-BFMI indicated no pathological behavior.

loo::pareto_k_table(loo(model_dist_apim_b))

All Pareto k estimates are good (k < 0.7).

Distinguishable Dyads - Check Model Convergence and Fit

plot(model_dist_apim_b, ask = FALSE)

Distinguishable Dyads - Check Model Convergence and Fit

pp_check(model_dist_apim_b, 'dens_overlay_grouped', group = 'is_male')

Distinguishable Dyads - Check Model Convergence and Fit

pp_check(model_dist_apim_b, 'ecdf_overlay_grouped', group = 'is_male')

Distinguishable Dyads - Check Model Convergence and Fit

pp_check(model_dist_apim_b, type = "loo_pit_overlay")

Distinguishable Dyads - Check Model Convergence and Fit

# Custom function to make DHARMa work with brms (see file 'Functions')
DHARMa.check_brms(model_dist_apim_b)

Distinguishable Dyads - Check Model Convergence and Fit

In case of non-convergence, bad residuals or low predictive accuracy, the model is likely miss-specified.

You can try different families that match your data by using the family() argument. For instance, you can easily conduct an ordinal regression by simply setting family = cumulative() and adjusting your priors. For dichotomous outcomes you can use family = bernoulli(). Brms supports a variety of families and link functions for these families. In general, you should use the default link function (check ?brms::family.brmsfit()).

Note that for some of the later models in this tutorial, some ESS values are too low. In this case, usually we would increase iterations (and warmup), and/or adapt_delta (see brms documentation).

Distinguishable Dyads - Results

Note that sigma is reported on the log scale (see ?brms::family.brmsfit()). Using exp(value) retrieves the simulated SD values well.

Distinguishable Dyads - (Intensive) Longitudinal APIM

Distinguishable Dyads - L-APIM - Simulated Dataset

Distinguishable Dyads - L-APIM - Centering

With repeated measures we need to center variables in order to not conflate levels of analysis:

• Between-Person: The person-mean in relation to the grand mean
• Within-Person: Daily fluctuations of an individual from their person-mean

While it may seem like we have 3-levels (couple/person/day), by including the means of both partners and their correlations in the model, all information about the couple-level (level 3) is included. This will become clear when comparing the exchangeable APIM to the DIM. We thus use a 2-level model with the dyad as the level of analysis.

Distinguishable Dyads - L-APIM - Centering

df_long_apim <- bmlm::isolate(
  df_long_apim, 
  by = 'userID', # NOT coupleID
  value = c('provided_support_actor', 'provided_support_partner'),
  which = 'both' 
) %>%
# renaming to avoid confusion with between- and within COUPLE centred DIM variables
  rename(
    provided_support_actor_cwp = provided_support_actor_cw,
    provided_support_partner_cwp = provided_support_partner_cw,
    provided_support_actor_cbp = provided_support_actor_cb,
    provided_support_partner_cbp = provided_support_partner_cb,
  ) %>%
  # Centering time
  mutate(
    diaryday_c = scale(diaryday, center = TRUE, scale = FALSE)
  )

Distinguishable Dyads - L-APIM - Model nlme

nlme model of this data for reference (simplified random effects structure to achieve convergence).

library(nlme)

model1 <- lme(
  fixed = 
    # Intercepts
    closeness ~ 0 + is_male + is_female +
    
    # Time Slopes
    is_male:diaryday_c +
    is_female:diaryday_c +
    
    # Between-Person APIM
    is_male:provided_support_actor_cbp + 
    is_male:provided_support_partner_cbp +
    
    is_female:provided_support_actor_cbp + 
    is_female:provided_support_partner_cbp +
  
    # Within-Person APIM
    is_male:provided_support_actor_cwp + 
    is_male:provided_support_partner_cwp +
    
    is_female:provided_support_actor_cwp + 
    is_female:provided_support_partner_cwp,
  
  random = 
    ~ 0 + is_male + is_female 
       
    + is_male:diaryday_c + is_female:diaryday_c
   
  # in nlme we need to remove many random slopes to achieve convergence
  # alternatively, we could try to remove correlations between the slopes. 
    #+ is_male:provided_support_actor_cwp 
    #+ is_male:provided_support_partner_cwp
    #+ is_female:provided_support_actor_cwp 
    #+ is_female:provided_support_partner_cwp 
  | coupleID, 
  
  weights = varIdent(form = ~ 1 | gender), # heterogeneous residual variances
  corr = corCompSymm(form = ~ 1 | coupleID/diaryday), # compound symmetry (partner's values on each day)
  
  data = df_long_apim, 
  control = list(maxIter=1000)
)

summary(model1)

Distinguishable Dyads - L-APIM - Model brms

Equivalent brms model with simplified random effects structure.

formula <- bf(
  closeness ~ 0 + is_male + is_female +
    
    # Time Slopes
    is_male:diaryday_c +
    is_female:diaryday_c +
    
    # Between-Person APIM
    is_male:provided_support_actor_cbp + 
    is_male:provided_support_partner_cbp +
    
    is_female:provided_support_actor_cbp + 
    is_female:provided_support_partner_cbp +
  
    # Within-Person APIM
    is_male:provided_support_actor_cwp + 
    is_male:provided_support_partner_cwp +
    
    is_female:provided_support_actor_cwp + 
    is_female:provided_support_partner_cwp +
  
  # Accounting for non-independence between partners' means and trajectories 
  # and effect sensitivities via random effects:
  (0 + is_male + is_female
    + is_male:diaryday_c + is_female:diaryday_c
  | coupleID) +
    
    
  # Accounting for daily non-independence (mimicking nlme's corCompSym)
  # modelled via common 'shocks' with a random intercept:
  (1 | coupleID:diaryday)

  # heteroscedastic residuals
  , sigma ~ 0 + is_male + is_female 
)

priors <- c(
  prior(normal(4, 2), class = "b", coef = "is_male"), # male intercept
  prior(normal(4, 2), class = "b", coef = "is_female"), # female intercept
  prior(normal(0, 2), class = "b") # other beta coefficients
  # Other priors can be specified, but brms choosed defaults well. 
  # Defaults can be seen via default_prior(formula, data)
)

model_dist_apim_long_simple <- brm(
  formula = formula, 
  data = df_long_apim,
  family = gaussian(link = identity), # student() is also often a nice robust option. 
  prior = priors,
  chains = 4,
  cores = 4,
  iter = 2000,
  warmup = 1000, 
  seed = 123,
  file = file.path('brms_cache', 'example1_dist_apim_long_simpel') # Cache the model
)

summary(model_dist_apim_long_simple)

Distinguishable Dyads - L-APIM - Results (scrollable)

Distinguishable Dyads - L-APIM - AR1 Model

We can try to add an AR1 residual structure.

    ... +
    ar(time = diaryday, gr = coupleID:userID, p = 1),
    
  sigma ~ 0 + is_male + is_female
)

Distinguishable Dyads - L-APIM - AR1 Model

Warning: Often, the daily shocks and the time slopes are enough and adding AR1 additionally may “soak up” too much variance. this can bias fixed effects estimates towards 0.

You may want to check which model performs better!

loo1 <- loo(model_dist_apim_long_simple)
loo2 <- loo(model_dist_apim_long_ar)

loo::pareto_k_table(loo1)

All Pareto k estimates are good (k < 0.7).

loo::pareto_k_table(loo2)

All Pareto k estimates are good (k < 0.7).

Distinguishable Dyads - L-APIM - AR1 Model

a <- loo_compare(loo1, loo2)
print(a)

                            elpd_diff se_diff
model_dist_apim_long_simple    0.0       0.0 
model_dist_apim_long_ar     -107.9      25.3 

report::report(a)

The difference in predictive accuracy, as indexed by Expected Log Predictive
Density (ELPD-LOO), suggests that 'model_dist_apim_long_simple' is the best
model (ELPD = -18689.18), followed by 'model_dist_apim_long_ar' (diff-ELPD =
-107.87 +- 25.27, p < .001)

In this instance, the AR1 model performs much worse, indicating potential overfitting or other issues. Even in the case of no clear difference, it may be advisable to choose the more parsimonious model, which is the one without AR1.

Maximal Model Specification (no AR(1))

We add all random slopes and their correlations. This is not converging in nlme.

formula <- bf(
  closeness ~ 
    
    # Intercepts
    0 + is_male + is_female +
    
    # Time Slopes
    is_male:diaryday_c +
    is_female:diaryday_c +
    
    # Between-Person APIM
    is_male:provided_support_actor_cbp + 
    is_male:provided_support_partner_cbp +
    
    is_female:provided_support_actor_cbp + 
    is_female:provided_support_partner_cbp +
  
    # Within-Person APIM
    is_male:provided_support_actor_cwp + 
    is_male:provided_support_partner_cwp +
    
    is_female:provided_support_actor_cwp + 
    is_female:provided_support_partner_cwp +
    
    # Accounting for non-independence between partners' means and trajectories 
    # and effect sensitivities via random effects:
    (0 + is_male + is_female +  
       
       is_male:diaryday_c + is_female:diaryday_c +
       
       is_male:provided_support_actor_cwp + 
       is_male:provided_support_partner_cwp +
       is_female:provided_support_actor_cwp + 
       is_female:provided_support_partner_cwp 
    | coupleID ) + 
    
    
    # Accounting for daily non-independence
    (1 | coupleID:diaryday), 
  
  # No more AR1. You could also test MA or ARMA. 
    
  sigma ~ 0 + is_male + is_female   # heteroscedastic residuals
)

priors <- c(
  prior(normal(4, 2), class = "b", coef = "is_male"),
  prior(normal(4, 2), class = "b", coef = "is_female"),
  prior(normal(0, 2), class = "b")
)


model_dist_apim_long_complex <- brm(
  formula = formula, 
  data = df_long_apim,
  family = gaussian(link = identity),
  prior = priors,
  chains = 4,
  cores = 4,
  iter = 2000,
  warmup = 1000, 
  seed = 123,
  file = file.path('brms_cache', 'example1_dist_apim_long_complex') # Cache the model
)

Maximal Model Specification: Results

Maximal Model Specification: Diagnostics

In complex models like this, we want to make sure we are not overfitting.

• Conduct regular model convergence checks (inspecting chains, Rhats, ESS, Multicollinearity, Residual diagnostics etc.). If these are bad (but good for a simpler model), it could be a sign of overfitting.

• Leave-one-out cross-validation is a powerful tool for investigating overfitting. It tests out of sample prediction and thus punishes overfitting.

Maximal Model Specification: Diagnostics

loo1 <- loo(model_dist_apim_long_simple)
loo2 <- loo(model_dist_apim_long_complex)

loo::pareto_k_table(loo1)

All Pareto k estimates are good (k < 0.7).

loo::pareto_k_table(loo2)

Pareto k diagnostic values:
                         Count Pct.    Min. ESS
(-Inf, 0.7]   (good)     4179  100.0%  299     
   (0.7, 1]   (bad)         1    0.0%  <NA>    
   (1, Inf)   (very bad)    0    0.0%  <NA>    

In case of High Pareto-Ks there are influential datapoints present. This can be a sign of overfitting but it is not necessarily a problem. But it is recommended to use loo(model, moment_match = TRUE) in these cases (the function will tell you).

The loo values (elpd) themselves cannot be interpreted individually, but we can compare models again (next page)

Maximal Model Specification: Diagnostics

a <- loo_compare(loo1, loo2)
print(a)

                             elpd_diff se_diff
model_dist_apim_long_simple   0.0       0.0   
model_dist_apim_long_complex -3.9       1.6   

No obvious overfitting detected. In fact, the more complex model seems to be more predictive. Some rule of thumbs say that if the difference is higher than 2x the SE or 4x the SE we can say there is a difference.

If there was no obvious difference, we could follow our theory and the philosophy of the maximal random effects structure (Barr et al., 2013). Alternatively, we could follow the philosophy of parsimony, especially in the territory of such complex models.

(continued on next page)

Maximal Model Specification: Diagnostics

One approach is to use the report package’s approach to obtain a p-value for the comparison (see their documentation for which assumptions they make):

report::report(a)

The difference in predictive accuracy, as indexed by Expected Log Predictive
Density (ELPD-LOO), suggests that 'model_dist_apim_long_simple' is the best
model (ELPD = -18689.18), followed by 'model_dist_apim_long_complex' (diff-ELPD
= -3.87 +- 1.59, p = 0.015)

In this instance, the maximal model (no AR) seems to be the preferred one.

In the case of overfitting, we could try to improve the model by subsequently excluding the smallest random effects (or correlations between them) and compare various models. For more guidance on random effects check out del Rosario & West (2025).

Maximal Model Specification: Comparing Estimates to Simulated Values (Ground Truth)

Exchangeable Dyads - Cross-Sectional APIM

(e.g., del Rosario & West, 2025; Kenny & Ackerman, 2023)

Exchangeable Dyads - Cross-Sectional APIM

(e.g., del Rosario & West, 2025; Kenny & Ackerman, 2023)

Exchangeable Dyads - Main Assumptions

Partners are exchangeable, i.e., not systematically different.

• Equal actor effects
• Equal partner effects
• Equal means
• Equal residual variances

But they should still be allowed to vary within each couple, while being correlated.

(e.g., del Rosario & West, 2025; Kenny & Ackerman, 2023)

Exchangeable Dyads - Cross-Sectional APIM: Data

Same cross-sectional data in the same 4-field actor-partner format as before.

print_df(head(df_apim))

Exchangeable Dyads - Cross-Sectional APIM: Model

formula <- bf(
  satisfaction ~ 1 + 
    communication_actor_gmc + communication_partner_gmc + 
    
    # Option 1: A single couple level random intercept.
    # In the cross sectional case this is the maximal random effects structure and
    # (1 | coupleID)
    
    # Option 2: using a compound symmetry residual structure. 
    cosy(gr = coupleID)
  
  # Note: no need to model separate sigmas for each partner.
  # Homogeneous residual variance is estimated. 
  # Implied: sigma = ~ 1
)

priors <- c(
  prior(normal(2, 10), class = "Intercept"),
  prior(normal(0, 5), class = "b"),
  prior(student_t(3, 0, 1.5), class = "sigma"),
  prior(beta(1, 3), class = "cosy")
)

model_ind_apim <- brm(
  formula = formula, 
  data = df_apim,
  family = gaussian(link = identity),
  prior = priors,
  chains = 4,
  cores = 4,
  iter = 2000,
  warmup = 1000, 
  seed = 123,
  file = file.path('brms_cache', 'example1_ind_apim') # Cache the model
)

Exchangeable Dyads - Cross-Sectional APIM: Results

Exchangeable Dyads - Cross-Sectional APIM: Test for Distinguishability

Leave-one-out (loo) cross-validation for model comparison (even if not nested).

a <- loo_compare(
  loo(model_ind_apim), 
  loo(model_dist_apim_b)
)
print(a)

report::report(a)

The difference in predictive accuracy, as indexed by Expected Log Predictive Density (ELPD-LOO), suggests that ‘model_dist_apim_b’ is the best model (ELPD = -1806.42), followed by ‘model_ind_apim’ (diff-ELPD = -172.85 +- 16.53, p < .001).
See: documentation of “report”

Exchangeable Dyads - Cross-Sectional DIM

Exchangeable Dyads - Cross-Sectional DIM: Data

Starting from scratch (no 4-field data needed)

Exchangeable Dyads - Cross-Sectional DIM: Centering

Decompose communication variance into:

• Between-couple (cbc): Couple-mean communication skills in relation to other couples
• Within-couple (cwc): Individuals’ communication skills in relation to their couple-mean

Same assumption about exchangeability as in the APIM

Exchangeable Dyads - Cross-Sectional DIM: Centering

df_dim <- df_apim2 %>%
  group_by(coupleID) %>%
  mutate(
    communication_cm = mean(communication, na.rm = TRUE),
    communication_cwc= communication - communication_cm
  ) %>%
  ungroup() %>%
  mutate(
    communication_cbc= communication_cm - mean(communication_cm, na.rm = TRUE)
  ) %>% 
  dplyr::select(c('userID', 'coupleID', 'satisfaction', 'communication_cwc', 'communication_cbc'))

print_df(head(df_dim))

Exchangeable Dyads - Cross-Sectional DIM: Model

formula <- bf(
  satisfaction ~ 1 + 
    communication_cbc+ communication_cwc + 
    cosy(gr = coupleID) 
)

priors <- c(
  prior(normal(2, 10), class = "Intercept"),
  prior(normal(0, 5), class = "b"),
  prior(student_t(3, 0, 1.5), class = "sigma"),
  prior(beta(1, 3), class = "cosy")
)

model_ind_dim <- brm(
  formula = formula, 
  data = df_dim,
  family = gaussian(link = identity),
  prior = priors,
  chains = 4,
  cores = 4,
  iter = 2000,
  warmup = 1000, 
  seed = 123,
  file = file.path('brms_cache', 'example1_ind_dim') # Cache the model
)

Exchangeable Dyads - Cross-Sectional DIM: Results

Equivalence APIM and DIM

\[b_{actor\_gmc} + b_{partner\_gmc} = b_{cbc}\] \[b_{actor\_gmc} - b_{partner\_gmc} = b_{cwc}\]

(Bolger et al., 2025)

Equivalence APIM and DIM

Equivalence APIM and DIM

• If the couple mean goes up by 1 and the within-couple deviation from the couple mean stays fixed, this means that both partners’ predictors must go up by 1 unit. Thus, the effects are linear combinations:

\[b_{cbc} = b_{actor\_gmc} + b_{partner\_gmc}\]

• Similarly, if a person’s deviation from their couple mean is one unit higher and the couple mean stays fixed, this means their partners’ deviation must be one unit lower. Thus, the effects are a linear combination:

\[b_{cwc} = b_{actor\_gmc} - b_{partner\_gmc}\]

(Bolger et al., 2025)

Equivalence APIM and DIM

• Conversely, to obtain the actor and partner effects given the DIM estimates:

\[b_{actor\_gmc} = \frac{b_{cbc} + b_{cwc}}{2}\]

\[b_{partner\_gmc} = \frac{b_{cbc} - b_{cwc}}{2}\]

Equivalence APIM and DIM

Example: using APIM coefficients to obtain DIM coefficients and computing Credible Intervals of DIM coefficients.

a <- hypothesis(
  model_ind_apim, 
  "communication_actor_gmc + communication_partner_gmc = 0"
)

round(a$hypothesis[,c(2,3,4,5)], 2)

  Estimate Est.Error CI.Lower CI.Upper
1     1.96      0.04     1.88     2.04

Equivalence APIM and DIM: Takeaway

• APIM and DIM are reparametrizations of the same model
• APIM: intuitive actor/partner framing
• DIM: clean separation of between vs within components
• Estimating one model allows for directly obtaining estimates of the other
• Same random-effects structure at dyad level and same assumptions

Side Note:

In the distinguishable case, the Dyadic Score Model (DSM) (e.g., Stadler et al., 2023) is equivalent to the APIM! (Bolger et al., 2025; Iida et al., 2018)

Exchangeable Dyads - Longitudinal APIM/DIM: Simulated Data

We use the same data as in the exchangeable case. But we don’t need the gender-related variables anymore, as we assume that partners are indistinguishable.

Exchangeable Dyads - Longitudinal DIM/APIM

• We need an APIM or DIM on both the between person level and the within-person level.
• Due to equivalence, we could have a DIM on one level and an APIM on the other.

If we want a within-person DIM:

df_long_dim <- df_long_apim2 %>%
  mutate(
    
    # Between Person Level DIM (centering between couples and within couple between person)
    provided_support_cbc = (provided_support_actor_cbp + provided_support_partner_cbp) / 2,
    provided_support_cwcbp = (provided_support_actor_cbp - provided_support_partner_cbp) / 2,
    
    # Within Person Level DIM (couple mean and deviation on each day)
    provided_support_cwp_mean = (provided_support_actor_cwp + provided_support_partner_cwp) / 2,
    provided_support_cwp_halfdiff = (provided_support_actor_cwp - provided_support_partner_cwp) / 2
  ) 

Exchangeable Dyads - Longitudinal DIM/APIM: Preparing the sum and difference approach

We need to randomly assign -1 and +1 for each Partner within each couple. This will be needed for contrast coding the intercept.

For example:

df_long_dim <- df_long_dim %>%
  group_by(coupleID) %>%
  mutate(
    base = ifelse(userID == min(userID), 1, -1),
    flip = 1 - 2 * rbinom(1, 1, 0.5),   # yields +1 or -1 once per couple
    Idiff = base * flip
  ) %>%
  ungroup() %>%
  relocate(Idiff, .after = coupleID) %>%
  dplyr::select(-base, -flip)


print_df(print_couple_preview(df_long_dim, couple_id = "31"))  

Exchangeable Dyads - Longitudinal DIM/APIM: Preparing the Sum and Difference Approach

Exchangeable Dyads - Longitudinal APIM/DIM: Model

formula <- bf(
  closeness ~ 1 + 
    
    diaryday_c +
  
    # Within-person APIM
    provided_support_actor_cwp + 
    provided_support_partner_cwp +
    # Equivalent to within-person DIM:
    # provided_support_cwp_mean + 
    # provided_support_cwp_halfdiff +
  
    # Between-person APIM
    provided_support_actor_cbp +
    provided_support_partner_cbp +
    # Equivalent to between-person DIM
    # provided_support_cwcbp+
    # provided_support_cbc +
    
    # Dyad-Level intercept and slopes for time-varying predictors
    (1 + diaryday_c + provided_support_actor_cwp + provided_support_partner_cwp | coupleID) +
    # Both partners' deviations from these dyad-level means and slopes.
    # Note: separate random effects block to make them uncorrelated form dyad-level REs. 
    (0 + Idiff + 
         I(Idiff * diaryday_c) + 
         I(Idiff * provided_support_actor_cwp) + 
         I(Idiff * provided_support_partner_cwp) | coupleID) +
    
    # Accounting for same-day shocks/coupling
    (1 | coupleID:diaryday)
    
    # Autocorrelated residuals
    # Due to bad fit in the distinguishable case, we omit AR1 from the start.
    # ar(time = diaryday, gr = coupleID:userID, p = 1)
  
  # Again, no need to model heterogeneous residual variances (sigma)
  # Implied: sigma = ~ 1
)

priors <- c(
  prior(normal(4, 2), class = "Intercept"),
  prior(normal(0, 2), class = "b")
  
  # We could set priors for other things, but brms sets good default priors
  #prior(student_t(3, 0, 1.5), class = "sd"),
  #prior(lkj(2), class = "cor"), # correlation prior for random effect matrix
  #prior(student_t(3, 0, 1.5), class = "sigma")
)

model_apim_ind_long <- brm(
  formula = formula, 
  data = df_long_dim,
  family = gaussian(link = identity),
  prior = priors,
  chains = 4,
  cores = 4,
  iter = 2000,
  warmup = 1000, 
  seed = 123,
  file = file.path('brms_cache', 'model_apim_ind_long_apim')
)

Exchangeable Dyads - Longitudinal APIM/DIM: Explanation Idiff

For appropriate random effects, we can use the Sum and Difference approach (del Rosario & West, 2025; Kenny & Ackerman, 2023):

• Randomly give one partner a constant of -1 and the other a constant of 1
• The couple level intercept represents the mean (or sum) of both partners’ intercepts. (1 | coupleID)
• A column of 1s and -1s represent deviations (difference) of each partner from the couple intercept with each partner contributing equally in one direction (+1 and -1)

Exchangeability is retained: if partners are flipped, the result is the same.

(del Rosario & West (2025) provide practical guidance on how to reduce the random effects structure in case of non-convergence.)

Exchangeable Dyads - Longitudinal APIM/DIM: Results

Comparing APIM (left) and DIM (right) output of fixed effects.

Equivalence holds on both levels. Numeric deviations arise due to priors, sampling noise, and rounding.

Exchangeable Dyads - Longitudinal APIM/DIM: Results

APIM Results in Detail:

Extract Full APIM Random Effects Variance-Covariance Matrix

We can rotate the random effects structure back to obtain the full APIM Random effects variance-covariance matrix (Kenny & Ackerman, 2023). Just as we would in a SEM model with equality constraints!

For example, the within-person variance of the intercept is:

\[var(Intercept) + var(Idiff)\]

and the cross-person covariance of both partners’ intercepts is:

\[var(Intercept) - var(Idiff)\]

Note: This only holds true when the means and Idiffs are uncorrelated. This is why we put them in separate random effects blocks.

Extract Full APIM Random Effects Variance-Covariance Matrix

# Custom function 
rotate_apim_covariance <- function(
    fit,
    gr = "coupleID",
    Idiff = "Idiff"
) {
  random_effects <- fit$ranef
  random_effects <- random_effects[random_effects$group == gr,]

  SUM <- random_effects$coef[!grepl(Idiff, random_effects$coef)]
  DIFF <- random_effects$coef[grepl(Idiff, random_effects$coef)]
  
  varcor <- VarCorr(fit)[[gr]]
  
  within_person_covariance_matrix <- varcor$cov[SUM,'Estimate',SUM] + varcor$cov[DIFF,'Estimate',DIFF]
  cross_person_covariance_matrix <- varcor$cov[SUM,'Estimate',SUM] - varcor$cov[DIFF,'Estimate',DIFF]
  
  p <- nrow(within_person_covariance_matrix)
  
  Full <- rbind(
    cbind(within_person_covariance_matrix,  cross_person_covariance_matrix),
    cbind(cross_person_covariance_matrix, within_person_covariance_matrix)
  )
  
  # nice labels
  base <- SUM
  rn <- c(paste0("PartnerA_", base), paste0("PartnerB_", base))
  colnames(Full) <- rn
  rownames(Full) <- rn
  
  return(list(within_person_covariance_matrix=within_person_covariance_matrix, cross_person_covariance_matrix=cross_person_covariance_matrix, full_covariance_matrix=Full))
}

a <- rotate_apim_covariance(model_apim_ind_long)

Extract Full APIM Random Effects Variance-Covariance Matrix: Within-Person Matrix

print_df(a$within_person_covariance_matrix)

Extract Full APIM Random Effects Variance-Covariance Matrix: Cross-Person Matrix

print_df(a$cross_person_covariance_matrix)

Extract Full APIM Random Effects Variance-Covariance Matrix: Full Matrix

print_df(a$full_covariance_matrix)

Convert to correlation matrix of random effects

# Custom function
sd_cor <- function(Sigma) {
  sds  <- sqrt(diag(Sigma))
  cors <- cov2cor(Sigma)  
  list(sd = sds, cor = cors)
}

full <- sd_cor(a$full_covariance_matrix)

Random effects SDs

print_df(as.data.frame(round(full$sd, 3)))

Random effect correlation matrix

print_df(as.data.frame(round(full$cor, 3)))